!Marche pas

!Dynamique moléculaire d'annihilation 1D avec condition initiale bimodale: \ln(N/N_now) \sim - \xi \ln(t) 
!But: transition entre régime en 1/\sqrt(t) [temps court - corrélations] et 1/t [temps longs - choas moléculaire]
!Version O(N)

module fonctions
implicit none
integer, parameter :: N = 6 !nombre initial de particules
double precision, parameter :: intervalle = 0.01d0 !intervalle log10(N/Nnow)
integer, parameter :: Navg = 200 !nombre de réalisations où toutes les particules sont annihilées pour moyenne
integer, parameter :: Nmin = 2 !nombre minimal de particules => arrête le run (laisser Nmin = 2)
character*(*), parameter :: filename = 'MDbimodal1d.txt' !nom du fichier de sortie
double precision, parameter :: pi = 3.1415926535897932384626433832795028841971693993751058209749445923078164
contains

FUNCTION ran2(idum)
INTEGER idum,IM1,IM2,IMM1,IA1,IA2,IQ1,IQ2,IR1,IR2,NTAB,NDIV
REAL :: ran2,AM,EPS,RNMX
PARAMETER (IM1=2147483563,IM2=2147483399,AM=1./IM1,IMM1=IM1-1,IA1=40014,IA2=40692,IQ1=53668,IQ2=52774,IR1=12211,&
          &IR2=3791,NTAB=32,NDIV=1+IMM1/NTAB,EPS=1.2e-7,RNMX=1.-EPS)
INTEGER :: idum2,j,k,iv(NTAB),iy
SAVE iv,iy,idum2
DATA idum2/123456789/, iv/NTAB*0/, iy/0/
if(idum.le.0)then
	idum=max(-idum,1)
	idum2=idum
	do j=NTAB+8,1,-1
		k=idum/IQ1
		idum=IA1*(idum-k*IQ1)-k*IR1
		if(idum.lt.0)idum=idum+IM1
		if(j.le.NTAB)iv(j)=idum
	enddo
	iy=iv(1)
endif
k=idum/IQ1
idum=IA1*(idum-k*IQ1)-k*IR1
if(idum.lt.0)idum=idum+IM1
k=idum2/IQ2
idum2=IA2*(idum2-k*IQ2)-k*IR2
if(idum2.lt.0)idum2=idum2+IM2
j=1+iy/NDIV
iy=iv(j)-idum2
iv(j)=idum
if(iy.lt.1)iy=iy+IMM1
ran2=min(AM*iy,RNMX)
return
end function ran2

SUBROUTINE Swop(A,B)
  DOUBLE PRECISION, DIMENSION(1:4), INTENT(INOUT) :: A, B
  DOUBLE PRECISION, DIMENSION(1:4)                :: Temp
  Temp = A
  A = B
  B = Temp
END SUBROUTINE Swop

RECURSIVE SUBROUTINE Quick_Sort(X,L,R)
  DOUBLE PRECISION, DIMENSION(1:N,1:4), INTENT(INOUT) :: X    ! list
  INTEGER, INTENT(IN)            		   								:: L, R ! left, right bounds
  INTEGER L1, R1
  IF (L < R)THEN
    L1 = L
    R1 = R
    DO
      DO WHILE (L1 < R .and. X(L1,1) <= X(L,1))  ! shift L1 right
        L1 = L1 + 1
      END DO
      DO WHILE (L < R1 .and. X(R1,1) >= X(L,1))  ! shift R1 left
        R1 = R1 - 1
      END DO
      IF (L1 < R1) CALL Swop(X(L1,:),X(R1,:)) ! swop
      IF (L1 >= R1) EXIT                     ! crossover - partition 
    END DO
    CALL Swop(X(L,:),X(R1,:))      ! partition with X(L) at R1          
    CALL Quick_Sort(X,L,R1-1) ! now attack left subproblem         
    CALL Quick_Sort(X,R1+1,R) ! don't forget right subproblem      
  END IF
END SUBROUTINE Quick_Sort

end module fonctions

!----------------------------------------------------------------------------------
!----------------------------------------------------------------------------------

program MDbimodal
use fonctions
implicit none
integer :: seed,i,k,avg,Nnow,taille,iintmin
double precision :: xold,dN,time,dNnow,tmp,tmp2,intmin
double precision, dimension(1:N,1:4) :: particules !particules(x,1) = position, particules(x,2) = vitesse
double precision, allocatable, dimension(:) :: x,t,nbmoyenne
seed = -14731
dN = dble(N)
taille = int(dlog10(dN/dble(Nmin))/intervalle) + 10 !première approximation
allocate(x(1:taille),t(1:taille),nbmoyenne(1:taille))
x=0.d0
t=0.d0
nbmoyenne = 0.d0
boucle_moyennes: do avg=1,Navg
	k = 1
	Nnow = N
	dNnow = dble(Nnow)
	time = 0.d0
	xold = 0.d0
	!***Initialisation: position initiale aléatoire \in [0,1]. Vitesse: +1 = vers la droite, -1 = vers la gauche
	particules = 0.d0
	do i=1,N
		particules(i,1) = dble(ran2(seed))
		if (ran2(seed) > 0.5d0) then
			particules(i,2) = 1.d0
		else
			particules(i,2) = -1.d0
		end if
	end do
	!***met dans l'ordre croissant les particules
	CALL Quick_Sort(particules,1,N)
	!***passe dans l'espace des distances relatives
	tmp = particules(1,1)
	tmp2 = particules(1,2)
	do i=1,N-1
		particules(i,1) = particules(i+1,1) - particules(i,1)
		particules(i,2) = particules(i+1,2)*particules(i,2)
	end do
	particules(N,1) = 1.d0-particules(N,1) + tmp
	particules(N,2) = particules(N,2)*tmp2
	do i=2,N-1
		particules(i,3) = dble(i-1)
		particules(i,4) = dble(i+1)
	end do
	particules(1,3) = dN
	particules(1,4) = 2.d0
	particules(N,3) = dble(N-1)
	particules(N,4) = 1.d0
	!***met dans l'ordre croissant les distances relatives
	CALL Quick_Sort(particules,1,N)

	boucle_un_run: do
		!***calcule la distance la plus courte (cond. bords périodiques)
		intmin = 10.d0

do i=1,N
write(*,*) particules(i,1),particules(i,2),nint(particules(i,3)),nint(particules(i,4))
end do
read(*,*)

		i=0
		pcoll: do
			if ((particules(i,2) < 0.d0) .and. (particules(i,1) <= intmin)) then
				exit pcoll
			else
				i = iintmin
				intmin = 10.d0
				CALL Quick_Sort(particules,1,N)
				exit pcoll
			end if
		end do pcoll	
		time = particules(i,1)

write(*,*) Nnow
read(*,*)
		
		!***annihile les particules et redéfinit les variables
		particules(nint(particules(i,3)),1) = particules(nint(particules(i,3)),1) + particules(i,1) + particules(nint(particules(i,4)),1)
		particules(nint(particules(i,3)),2) = - particules(nint(particules(i,3)),2)*particules(nint(particules(i,4)),2)
		if (nint(particules(i,4)) /= N) then
			particules(nint(particules(i,3)),4) = particules(i,4) + 1.d0
		else
			particules(nint(particules(i,3)),4) = 1.d0
		end if
		if (nint(particules(i,4)) /= N) then
			particules(nint(particules(i,4)+1),3) = particules(i,3)
		else
			particules(1,3) = particules(i,3)
		end if
		particules(i,1) = 10.d0
		Nnow = Nnow - 2
		dNnow = dble(Nnow)
		!***enregistre la valeur minimale du nouvel intervalle créé
		if (particules(nint(particules(i,3)),1) < intmin) then
			intmin = particules(particules(i,3),1)
			iintmin = i
			CALL Quick_Sort(particules,1,N)
		end if
		!***répartition logarithmique des valeurs
		if (dlog10(dN/dNnow)-xold >= intervalle) then
			xold = dlog10(dN/dNnow)
			k = k+1
			x(k) = x(k) + xold
			t(k) = t(k) + time
			nbmoyenne(k)=nbmoyenne(k)+1
!			open(unit=1,file="state")
			write(*,'(" Run:",i5,"/",i5,"   Nnow: ",i7)') avg,Navg,Nnow
!			close(1)
		end if
		!***fin du run: écriture de toutes les valeurs, premier terme t=0 à éliminer car on a \ln(t) = \ln(0)
		if (Nnow <= Nmin) then
			if ((k > taille) .or. (avg==1)) taille = k	
			open(unit=1, file=filename)
			do i=2,taille
				write(1,*) x(i)/dble(nbmoyenne(i)),dlog10(t(i)/(dble(nbmoyenne(i))))
			end do
			close(1)
			exit boucle_un_run
		end if
	end do boucle_un_run

end do boucle_moyennes
deallocate(x,t,nbmoyenne)
write(*,*)
write(*,*) "Programme termine."
read(*,*)
end program MDbimodal